1 Introducción


1.1 Presentación

Esta práctica cubre de forma transversal la asignatura.

Las Prácticas 1 y 2 de la asignatura se plantean de una forma conjunta de modo que la Práctica 2 será continuación de la 1.

El objetivo global de las dos prácticas consiste en seleccionar uno o varios juegos de datos, realizar las tareas de preparación y análisis exploratorio con el objetivo de disponer de datos listos para aplicar algoritmos de clustering, asociación y clasificación.

1.2 Competencias

Las competencias que se trabajan en esta prueba son:

  • Uso y aplicación de las TIC en el ámbito académico y profesional.
  • Capacidad para innovar y generar nuevas ideas.
  • Capacidad para evaluar soluciones tecnológicas y elaborar propuestas de proyectos teniendo en cuenta los recursos, las alternativas disponibles y las condiciones de mercado.
  • Conocer las tecnolog??as de comunicaciones actuales y emergentes as?? como saberlas aplicar convenientemente para diseñar y desarrollar soluciones basadas en sistemas y tecnolog??as de la información.
  • Aplicación de las técnicas espec??ficas de ingenier??a del software en las diferentes etapas del ciclo de vida de un proyecto.
  • Capacidad para aplicar las técnicas espec??ficas de tratamiento, almacenamiento y administración de datos.
  • Capacidad para proponer y evaluar diferentes alternativas tecnológicas para resolver un problema concreto.

1.3 Objetivos

La correcta asimilación de todos los aspectos trabajados durante el semestre.
En esta práctica abordamos un caso real de miner??a de datos donde tenemos que poner en juego todos los conceptos trabajados. Hay que trabajar todo el ciclo de vida del proyecto. Desde el objetivo del proyecto hasta la implementación del conocimiento encontrado pasando por la preparación, limpieza de los datos, conocimiento de los datos, generación del modelo, interpretación y evaluación.

1.4 Descripción de la PEC a realizar

1.5 Recursos Básicos

Material docente proporcionado por la UOC.

1.6 Criterios de valoración

Ejercicios prácticos

Para todas las PEC es necesario documentar en cada apartado del ejercicio práctico que se ha hecho y como se ha hecho.

1.7 Formato y fecha de entrega

El formato de entrega es: usernameestudiante-PECn.html/doc/docx/odt/pdf
Fecha de entrega: 15/01/2020
Se debe entregar la PEC en el buzón de entregas del aula

1.8 Nota: Propiedad intelectual

A menudo es inevitable, al producir una obra multimedia, hacer uso de recursos creados por terceras personas. Es por lo tanto comprensible hacerlo en el marco de una práctica de los estudios de Informática, Multimedia y Telecomunicación de la UOC, siempre y cuando esto se documente claramente y no suponga plagio en la práctica.

Por lo tanto, al presentar una práctica que haga uso de recursos ajenos, se debe presentar junto con ella un documento en que se detallen todos ellos, especificando el nombre de cada recurso, su autor, el lugar donde se obtuvo y su estatus legal: si la obra esta protegida por el copyright o se acoge a alguna otra licencia de uso (Creative Commons, licencia GNU, GPL …). El estudiante deberá asegurarse de que la licencia no impide espec??ficamente su uso en el marco de la práctica. En caso de no encontrar la información correspondiente tendrá que asumir que la obra esta protegida por copyright.

Deberéis, además, adjuntar los ficheros originales cuando las obras utilizadas sean digitales, y su código fuente si corresponde.


2 Enunciado


Como continuación del estudio iniciado en la práctica 1, procedemos en esta práctica 2 a aplicar modelos anal??ticos sobe el juego de datos seleccionado y preparado en la práctica anterior.

De este modo se pide al estudiante que complete los siguientes pasos:

  1. Aplicar un modelo de generación de reglas a partir de Reglas de asociacion (Hay que usar el algoritmo arules).

  2. Aplicar un modelo no supervisado y basado en el concepto de distancia, sobre el juego de datos (Hay que usar el algoritmo kmeans).

  3. Aplica de nuevo el modelo anterior, pero usando una métrica distinta y compara los resultados (Hay que usar un algoritmo de clustering diferente al kmeans)

  4. Aplicar un modelo supervisado sobre el juego de datos sin haber aplicado previamente PCA/SVD (Hay que usar un árbol de decisión).

  5. Aplicar un modelo supervisado sobre el juego de datos habiendo aplicado previamente PCA/SVD (Hay que usar un árbol de decisión tras la aplicación de un PCA/SVD).

  6. ¿Ha habido mejora en capacidad predictiva, tras aplicar PCA/SVD? ¿A qué crees que es debido? (Hay que comparar los resultados de los puntos 4 y 5).


3 Solución


3.1 Carga de los datos

Vamos a comenzar realizando la carga de las librerías y funciones necesarias para ejecutar el ejercicio.

libraries <- c("caret", "arules", "cluster", "fpc","tibble", "factoextra", "C50", "ggbiplot", "dplyr", "cluster")

check.libraries <- is.element(libraries, installed.packages()[, 1]) == FALSE
libraries.to.install <- libraries[check.libraries]

if (length(libraries.to.install != 0)) {
  install.packages(libraries.to.install)
  }

success <- sapply(libraries, require, quietly = FALSE, character.only = TRUE)
if(length(success) != length(libraries)) {
  stop("A package failed to return a success in require() function.")
  }
normalize <- function(x) {
  return ((x - min(x)) / (max(x) - min(x)))
}

De la resolución de la práctica anterior, se han exportado los datos finales y vueltos a importar dentro del contexto de este proyecto.

data.raw <- read.csv('dengue_data.csv',stringsAsFactors = FALSE)
data.raw$cases <- as.factor(data.raw$cases)
summary(data.raw)
##        X                ID             city                year     
##  Min.   :   1.0   Min.   :   1.0   Length:1456        Min.   :1990  
##  1st Qu.: 364.8   1st Qu.: 364.8   Class :character   1st Qu.:1997  
##  Median : 728.5   Median : 728.5   Mode  :character   Median :2002  
##  Mean   : 728.5   Mean   : 728.5                      Mean   :2001  
##  3rd Qu.:1092.2   3rd Qu.:1092.2                      3rd Qu.:2005  
##  Max.   :1456.0   Max.   :1456.0                      Max.   :2010  
##    weekofyear     trimester              ndvi        precipitation   
##  Min.   : 1.00   Length:1456        Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:13.75   Class :character   1st Qu.:0.3239   1st Qu.:0.0255  
##  Median :26.50   Mode  :character   Median :0.4085   Median :0.0991  
##  Mean   :26.50                      Mean   :0.4356   Mean   :0.1172  
##  3rd Qu.:39.25                      3rd Qu.:0.5353   3rd Qu.:0.1793  
##  Max.   :53.00                      Max.   :1.0000   Max.   :1.0000  
##  air_temp_mean    diur_temp_range      humidity       total_cases    
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :  0.00  
##  1st Qu.:0.4005   1st Qu.:0.06621   1st Qu.:0.4397   1st Qu.:  5.00  
##  Median :0.5323   Median :0.10273   Median :0.6116   Median : 12.00  
##  Mean   :0.5375   Mean   :0.24174   Mean   :0.5752   Mean   : 24.68  
##  3rd Qu.:0.6864   3rd Qu.:0.42235   3rd Qu.:0.7157   3rd Qu.: 28.00  
##  Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :461.00  
##       cases     
##  low     :1431  
##  moderate:  15  
##  severe  :  10  
##                 
##                 
## 
str(data.raw)
## 'data.frame':    1456 obs. of  13 variables:
##  $ X              : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ ID             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ city           : chr  "sj" "sj" "sj" "sj" ...
##  $ year           : int  1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
##  $ weekofyear     : int  18 19 20 21 22 23 24 25 26 27 ...
##  $ trimester      : chr  "T2" "T2" "T2" "T2" ...
##  $ ndvi           : num  0.408 0.419 0.379 0.506 0.556 ...
##  $ precipitation  : num  0.0318 0.0584 0.0884 0.0393 0.0193 ...
##  $ air_temp_mean  : num  0.388 0.473 0.548 0.575 0.646 ...
##  $ diur_temp_range: num  0.0867 0.0691 0.0643 0.073 0.113 ...
##  $ humidity       : num  0.263 0.418 0.587 0.567 0.628 ...
##  $ total_cases    : int  4 5 4 3 6 2 4 5 10 6 ...
##  $ cases          : Factor w/ 3 levels "low","moderate",..: 1 1 1 1 1 1 1 1 1 1 ...

Dado que estos datos provienen de un ejercicio en el que se ha elaborado una amplia y exhaustiva tarea de limpieza y normalización de datos, vamos a saltarnos este paso y pasar directamente a la elaboración de los modelos.

Sin embargo, existe un problema que debe ser solventado y que ha sido descubierto durante la fase preliminar, y no documentada, de la elaboración del ejercicio. El problema surgió tras explorar la generación de reglas mediante arules cuando el algoritmo no generaba regla alguna, o lo que es lo mismo, todas las reglas estaban orientadas a la misma clase resultado.

Tras múltiples intentos y evaluaciones se observó que la variable resultante cases tiene una distribución poco equilibrada y no apta para la generación de reglas. Como era de esperar, este problema no estaba presente en los modelos no supervisados, aunque todos los modelos supervisamos daban como resultado una precisión superior al 98%, algo poco común.

Aunque tiene sentido agrupar los casos de dengue de forma que a priori nos parezca razonable, la realidad es que desde un punto de vista estadístico no lo es. Si observamos, más del 95% de los casos de dengue caen bajo la categoría de low o bajos. Esto hace que cualquier modelo generado a partir de estos datos sea poco o nada eficaz. El modelo será muy eficiente prediciendo la baja indicencia de dengue pero nefasto en la predicción de casos moderados o de alta incidencia de la enfermedad.

Para resolver este problema de distribución de las clases resultantes, se ha optado por la redistribución de esta variable siguiendo otros criterios más objetivos. Para ello se ha hecho uso de la media como línea de corte; es decir, que todas aquellas semanas donde se hayan detectado hasta 25 casos de dengue se considerará una semana normal, mientras que si superamos los 25 casos lo catalogaremos como pandemic.

summary(data.raw$total_cases)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   12.00   24.68   28.00  461.00

Tras realizar la nueva discretización de los datos mediante el método fixed, objetemos una distribución de 1036 casos normales y 420 pandémicos. Esta nueva distribución puede, sin duda alguna, generar modelos más fiables en su predicción.

outcome <- data.frame('outcome' = discretize(data.raw$total_cases,method="fixed", breaks = c(0,25, 500), labels=c('normal', 'pandemic')))
data.raw <- add_column(data.raw, outcome, .after = 13)
data.raw$cases <- NULL
summary(data.raw)
##        X                ID             city                year     
##  Min.   :   1.0   Min.   :   1.0   Length:1456        Min.   :1990  
##  1st Qu.: 364.8   1st Qu.: 364.8   Class :character   1st Qu.:1997  
##  Median : 728.5   Median : 728.5   Mode  :character   Median :2002  
##  Mean   : 728.5   Mean   : 728.5                      Mean   :2001  
##  3rd Qu.:1092.2   3rd Qu.:1092.2                      3rd Qu.:2005  
##  Max.   :1456.0   Max.   :1456.0                      Max.   :2010  
##    weekofyear     trimester              ndvi        precipitation   
##  Min.   : 1.00   Length:1456        Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:13.75   Class :character   1st Qu.:0.3239   1st Qu.:0.0255  
##  Median :26.50   Mode  :character   Median :0.4085   Median :0.0991  
##  Mean   :26.50                      Mean   :0.4356   Mean   :0.1172  
##  3rd Qu.:39.25                      3rd Qu.:0.5353   3rd Qu.:0.1793  
##  Max.   :53.00                      Max.   :1.0000   Max.   :1.0000  
##  air_temp_mean    diur_temp_range      humidity       total_cases    
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :  0.00  
##  1st Qu.:0.4005   1st Qu.:0.06621   1st Qu.:0.4397   1st Qu.:  5.00  
##  Median :0.5323   Median :0.10273   Median :0.6116   Median : 12.00  
##  Mean   :0.5375   Mean   :0.24174   Mean   :0.5752   Mean   : 24.68  
##  3rd Qu.:0.6864   3rd Qu.:0.42235   3rd Qu.:0.7157   3rd Qu.: 28.00  
##  Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :461.00  
##      outcome    
##  normal  :1036  
##  pandemic: 420  
##                 
##                 
##                 
## 

Normalizamos la columna weekofyears para poder hacer uso de esta variable en un futuro como variable discreata.

weekofyear_norm <-  as.data.frame(apply(data.raw[which( colnames(data.raw)=="weekofyear" )], 2, normalize))
data.raw$weekofyear <- NULL
data.raw <- add_column(data.raw, weekofyear_norm, .after = 4)
summary(data.raw)
##        X                ID             city                year     
##  Min.   :   1.0   Min.   :   1.0   Length:1456        Min.   :1990  
##  1st Qu.: 364.8   1st Qu.: 364.8   Class :character   1st Qu.:1997  
##  Median : 728.5   Median : 728.5   Mode  :character   Median :2002  
##  Mean   : 728.5   Mean   : 728.5                      Mean   :2001  
##  3rd Qu.:1092.2   3rd Qu.:1092.2                      3rd Qu.:2005  
##  Max.   :1456.0   Max.   :1456.0                      Max.   :2010  
##    weekofyear      trimester              ndvi        precipitation   
##  Min.   :0.0000   Length:1456        Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.2452   Class :character   1st Qu.:0.3239   1st Qu.:0.0255  
##  Median :0.4904   Mode  :character   Median :0.4085   Median :0.0991  
##  Mean   :0.4905                      Mean   :0.4356   Mean   :0.1172  
##  3rd Qu.:0.7356                      3rd Qu.:0.5353   3rd Qu.:0.1793  
##  Max.   :1.0000                      Max.   :1.0000   Max.   :1.0000  
##  air_temp_mean    diur_temp_range      humidity       total_cases    
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :  0.00  
##  1st Qu.:0.4005   1st Qu.:0.06621   1st Qu.:0.4397   1st Qu.:  5.00  
##  Median :0.5323   Median :0.10273   Median :0.6116   Median : 12.00  
##  Mean   :0.5375   Mean   :0.24174   Mean   :0.5752   Mean   : 24.68  
##  3rd Qu.:0.6864   3rd Qu.:0.42235   3rd Qu.:0.7157   3rd Qu.: 28.00  
##  Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :461.00  
##      outcome    
##  normal  :1036  
##  pandemic: 420  
##                 
##                 
##                 
## 

Como paso extra, vamos a buscar correlaciones en los datos para evaluar si es posible deshacernos de alguna variable. Los datos indican que no existe correlación entre los datos ya que ninguno es próximo a -1 o 1.

cor(data.raw[,c(5,7,8,9,10,11)])
##                 weekofyear        ndvi precipitation air_temp_mean
## weekofyear      1.00000000  0.04094926    0.11662773    0.42134406
## ndvi            0.04094926  1.00000000    0.17711009   -0.30325048
## precipitation   0.11662773  0.17711009    1.00000000   -0.01721905
## air_temp_mean   0.42134406 -0.30325048   -0.01721905    1.00000000
## diur_temp_range 0.07712561  0.67396605    0.20132179   -0.28037539
## humidity        0.34718816  0.07904915    0.45060559    0.50838016
##                 diur_temp_range   humidity
## weekofyear           0.07712561 0.34718816
## ndvi                 0.67396605 0.07904915
## precipitation        0.20132179 0.45060559
## air_temp_mean       -0.28037539 0.50838016
## diur_temp_range      1.00000000 0.01294761
## humidity             0.01294761 1.00000000

3.2 Generación de reglas

Para generar las reglas vamos a trabajar con dos subconjuntos de datos. En el primero utilizaremos la variable trimestre, que como ya sabemos, es una indicadora de las estaciones del año, y en el segundo trabajaremos sólo con datos climatológicos. La razón por la que hemos decidido utilizar dos subconjuntos de datos es porque creemos que el dato temporal es, en cierto modo, una función de los datos climatológicos y morfológicos recogidos.

Queremos ver la influencia directa que la variable trimestre pueda tener sobre las reglas en contraste con los factores climatológicos y morfológicos por si solos, y lo más importante, estudiar la eficacia de los modelos generados de ambos conjuntos.

3.2.1 Uso de variables estacionales

data.arules.seasons <- data.raw[, -c(1, 2, 3, 4, 5, 12)]
summary(data.arules.seasons)
##   trimester              ndvi        precipitation    air_temp_mean   
##  Length:1456        Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  Class :character   1st Qu.:0.3239   1st Qu.:0.0255   1st Qu.:0.4005  
##  Mode  :character   Median :0.4085   Median :0.0991   Median :0.5323  
##                     Mean   :0.4356   Mean   :0.1172   Mean   :0.5375  
##                     3rd Qu.:0.5353   3rd Qu.:0.1793   3rd Qu.:0.6864  
##                     Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##  diur_temp_range      humidity          outcome    
##  Min.   :0.00000   Min.   :0.0000   normal  :1036  
##  1st Qu.:0.06621   1st Qu.:0.4397   pandemic: 420  
##  Median :0.10273   Median :0.6116                  
##  Mean   :0.24174   Mean   :0.5752                  
##  3rd Qu.:0.42235   3rd Qu.:0.7157                  
##  Max.   :1.00000   Max.   :1.0000

A continuación, discretizamos las variables numéricas.

data.arules.seasons$ndvi <- discretize(data.arules.seasons$ndvi, method="interval", labels=c('water', 'cement', 'vegetal'))
data.arules.seasons$precipitation <- discretize(data.arules.seasons$precipitation, method="interval", labels=c('dry', 'normal', 'storm'))
data.arules.seasons$air_temp_mean <- discretize(data.arules.seasons$air_temp_mean, method="interval", labels=c('low', 'medium', 'high'))
data.arules.seasons$diur_temp_range <- discretize(data.arules.seasons$diur_temp_range, method="interval", labels=c('low', 'medium', 'high'))
data.arules.seasons$humidity <- discretize(data.arules.seasons$humidity, method="interval", labels=c('dry', 'normal', 'wet'))
summary(data.arules.seasons)
##   trimester              ndvi     precipitation air_temp_mean diur_temp_range
##  Length:1456        water  :399   dry   :1398   low   :184    low   :1000    
##  Class :character   cement :919   normal:  55   medium:872    medium: 332    
##  Mode  :character   vegetal:138   storm :   3   high  :400    high  : 124    
##    humidity       outcome    
##  dry   :158   normal  :1036  
##  normal:738   pandemic: 420  
##  wet   :560

Creamos los conjuntos de entramiento y test.

set.seed(777)
train <- createDataPartition(data.arules.seasons$outcome, p = 0.66, list = FALSE)

data.arules.seasons.train <- data.arules.seasons[train,]
data.arules.seasons.test <- data.arules.seasons[-train,]

data.arules.seasons.train.x <- data.arules.seasons.train[,1:6]
data.arules.seasons.train.y <- data.arules.seasons.train[,7]

data.arules.seasons.test.x <- data.arules.seasons.test[,1:6]
data.arules.seasons.test.y <- data.arules.seasons.test[,7]

Y estudiamos los resultados del modelo generado.

data.arules.seasons.model <- C5.0(x = data.arules.seasons.train.x, y = data.arules.seasons.train.y, rules=TRUE)
summary(data.arules.seasons.model)
## 
## Call:
## C5.0.default(x = data.arules.seasons.train.x, y =
##  data.arules.seasons.train.y, rules = TRUE)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Mon Jun  8 20:53:17 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 962 cases (7 attributes) from undefined.data
## 
## Rules:
## 
## Rule 1: (290/19, lift 1.3)
##  diur_temp_range in {medium, high}
##  ->  class normal  [0.932]
## 
## Rule 2: (482/78, lift 1.2)
##  trimester in {T2, T1}
##  ->  class normal  [0.837]
## 
## Rule 3: (322/132, lift 2.0)
##  trimester in {T3, T4}
##  diur_temp_range = low
##  ->  class pandemic  [0.590]
## 
## Default class: normal
## 
## 
## Evaluation on training data (962 cases):
## 
##          Rules     
##    ----------------
##      No      Errors
## 
##       3  220(22.9%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##     552   132    (a): class normal
##      88   190    (b): class pandemic
## 
## 
##  Attribute usage:
## 
##   83.58% trimester
##   63.62% diur_temp_range
## 
## 
## Time: 0.0 secs

Las reglas encontradas indican que:

Rango en temperatura diurna (medio, alto) -> casos normales

Estaciones de invierno y primavera -> casos normales

Estaciones verano y otoño + Rango en temperatura diurna (bajo) -> casos pandémicos.

Según estas reglas, existen dos factores que influencian la proliferación de los casos de dengue. Estos son las estaciones del año y la temperatura media diurna. Los casos aumentan durante el verano y el otoño (estación de huracanes para el área de San juan) y con un rango de temperatura media bajo; esto es, con pocas o ninguna bajada de temperatura durante la tarde. Los casos normales de dengue ocurren bajo condiciones totalmente opuestas. Podemos decir que al mosquito le gusta el calor y el agua.

A pesar de ello, no vemos que el modelo pueda describir muy bien las condiciones exactas bajo las que prolifera el mosquito ya que sólo describe el 50% de los casos pandémicos. También podemos observar que el modelo deja la predicción en manos de la estación del año ignorando el resto de las variables.

data.arules.seasons.model.chart <- C5.0(x = data.arules.seasons.train.x, y = data.arules.seasons.train.y)
plot(data.arules.seasons.model.chart, subtree = 1)

3.2.2 Uso de variables no estacionales

Vamos a proceder ahora con el segundo conjunto de datos donde ignoramos la estación del año.

data.arules.no_seasons <- data.arules.seasons[,-1]
summary(data.arules.no_seasons)
##       ndvi     precipitation air_temp_mean diur_temp_range   humidity  
##  water  :399   dry   :1398   low   :184    low   :1000     dry   :158  
##  cement :919   normal:  55   medium:872    medium: 332     normal:738  
##  vegetal:138   storm :   3   high  :400    high  : 124     wet   :560  
##      outcome    
##  normal  :1036  
##  pandemic: 420  
## 

Creamos los conjuntos de entramiento y test.

set.seed(777)
train <- createDataPartition(data.arules.no_seasons$outcome, p = 0.66, list = FALSE)

data.arules.no_seasons.train <- data.arules.no_seasons[train,]
data.arules.no_seasons.test <- data.arules.no_seasons[-train,]

data.arules.no_seasons.train.x <- data.arules.no_seasons.train[,1:5]
data.arules.no_seasons.train.y <- data.arules.no_seasons.train[,6]

data.arules.no_seasons.test.x <- data.arules.no_seasons.test[,1:5]
data.arules.no_seasons.test.y <- data.arules.no_seasons.test[,6]

Y estudiamos los resultados del modelo generado.

data.arules.no_seasons.model <- C5.0(x = data.arules.no_seasons.train.x, y = data.arules.no_seasons.train.y, rules=TRUE)
summary(data.arules.no_seasons.model)
## 
## Call:
## C5.0.default(x = data.arules.no_seasons.train.x, y
##  = data.arules.no_seasons.train.y, rules = TRUE)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Mon Jun  8 20:53:17 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 962 cases (6 attributes) from undefined.data
## 
## Rules:
## 
## Rule 1: (290/19, lift 1.3)
##  diur_temp_range in {medium, high}
##  ->  class normal  [0.932]
## 
## Rule 2: (693/147, lift 1.1)
##  air_temp_mean in {low, medium}
##  ->  class normal  [0.787]
## 
## Rule 3: (254/124, lift 1.8)
##  air_temp_mean = high
##  diur_temp_range = low
##  ->  class pandemic  [0.512]
## 
## Default class: normal
## 
## 
## Evaluation on training data (962 cases):
## 
##          Rules     
##    ----------------
##      No      Errors
## 
##       3  272(28.3%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##     560   124    (a): class normal
##     148   130    (b): class pandemic
## 
## 
##  Attribute usage:
## 
##   98.44% air_temp_mean
##   56.55% diur_temp_range
## 
## 
## Time: 0.0 secs

Aquí obtenemos las siguientes reglas:

Rango en temperatura diurna (medio, alto) -> casos normales

Temperatura aire promedio (bajo, medio) -> casos normales

Rango en temperatura diurna (bajo) + Temperatura aire promedio (alto) -> casos pandémicos

Para este conjunto de datos, el árbol resultante es más diverso al proveer reglas que incluyen otras variables, podemos decir que “se explica mejor” que el modelo anterior. Aquí podemos ver que los casos de dengue aumentan con las combinaciones de altas temperaturas promedios diarias y constantes durante la tarde y sobre todo en zonas urbanas (ndvi ~ 0).

El factor temperatura ha sido constante en ambos estudios, extendido los resultados a la temperatura media en este segundo conjunto de datos donde además se ha sumado el factor ndvi en forma de zonas con cemento o urbanizadas, que es precisamente donde el agua es más propensa a estancarse. A esta informacion, y según el modelo anterior, podemos añadirle el hecho de que las estaciones de verano y otoño son las mas propensas para desarrollo del mosquito.

data.arules.no_seasons.model.chart <- C5.0(x = data.arules.no_seasons.train.x, y = data.arules.no_seasons.train.y)
plot(data.arules.no_seasons.model.chart, subtree = 1)

Vamos a continuación a comparar ambos modelos en términos de eficiencia predictiva de los datos de prueba.

data.arules.seasons.predicted_model <- predict(data.arules.seasons.model, data.arules.seasons.test.x, type="class")
print(sprintf("La precisión del árbol con estaciones es: %.4f %%",100*sum(data.arules.seasons.predicted_model == data.arules.seasons.test.y) / length(data.arules.seasons.predicted_model)))
## [1] "La precisión del árbol con estaciones es: 78.9474 %"
data.arules.no_seasons.predicted_model <- predict(data.arules.no_seasons.model, data.arules.no_seasons.test.x, type="class")
print(sprintf("La precisión del árbol sin estaciones es: %.4f %%",100*sum(data.arules.no_seasons.predicted_model == data.arules.no_seasons.test.y) / length(data.arules.no_seasons.predicted_model)))
## [1] "La precisión del árbol sin estaciones es: 71.2551 %"

El árbol deducido del conjunto de datos que incluye las estaciones del año ha devengado un mayor porcentaje en la precisión de los resultados en contraste con aquel que utiliza datos puramente climáticos y morfológicos. Por el contrario, la definición de reglas es más rica y explicita al omitir la variable estacional.

3.3 Modelo no supervisado

Para realizar el modelo no supervisado vamos a repetir el procedimiento anterior y hacer uso de un conjunto con variable temporal y otro sin ella. Utilizaremos la variable normalizada weekofyear_morm. El propósito de este ejercicio es encontrar el número de clústeres óptimo para nuestro modelo mediante el método kmeans.

3.3.1 Uso de variables estacionales

data.kmeans.seasons <- data.raw[, c(5,7,8,9,10,11,13)]
summary(data.kmeans.seasons)
##    weekofyear          ndvi        precipitation    air_temp_mean   
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.2452   1st Qu.:0.3239   1st Qu.:0.0255   1st Qu.:0.4005  
##  Median :0.4904   Median :0.4085   Median :0.0991   Median :0.5323  
##  Mean   :0.4905   Mean   :0.4356   Mean   :0.1172   Mean   :0.5375  
##  3rd Qu.:0.7356   3rd Qu.:0.5353   3rd Qu.:0.1793   3rd Qu.:0.6864  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##  diur_temp_range      humidity          outcome    
##  Min.   :0.00000   Min.   :0.0000   normal  :1036  
##  1st Qu.:0.06621   1st Qu.:0.4397   pandemic: 420  
##  Median :0.10273   Median :0.6116                  
##  Mean   :0.24174   Mean   :0.5752                  
##  3rd Qu.:0.42235   3rd Qu.:0.7157                  
##  Max.   :1.00000   Max.   :1.0000

Para el cálculo de clusters óptimo, realizamos una iteración en la que asumimos un total de 2 hasta de 10 clusters. Con cada uno de estos valores utilizamos el algoritmo kmeans para clasificar las muestras y el método silhouette para estudiar la precisión en la agrupación de los datos.

x <- data.kmeans.seasons[,1:6]
d <- daisy(x) 

resultados <- rep(0, 7)
cluster  <- data.frame(c(2,3,4,5,6,7))

for (i in c(2,3,4,5,6,7))
{
  fit           <- kmeans(x, i)
  y_cluster     <- fit$cluster
  sk            <- silhouette(y_cluster, d)
  resultados[i] <- mean(sk[,3])
}

plot(2:7,resultados[2:7],type="o",col="blue",pch=0,xlab="Número de clusters",ylab="Silueta")

La gráfica revela una agrupación más definida para un modelo con 4 clústeres. El método elbow, debajo, recomienda el uso de 3 clústeres.

resultados <- rep(0, 7)
for (i in c(2,3,4,5,6,7))
{
  fit           <- kmeans(x, i)
  resultados[i] <- fit$tot.withinss
}
plot(2:7,resultados[2:7],type="o",col="blue",pch=0,xlab="Número de clusters",ylab="tot.tot.withinss")

Si procedemos con la función kmeansruns bajos los criterios de la silueta media y Calinski-Harabasz obtenemos que para ambos casos el número óptimo de clusters es de 4.

fit_ch  <- kmeansruns(x, krange = 1:10, criterion = "ch") 
fit_asw <- kmeansruns(x, krange = 1:10, criterion = "asw")
fit_ch$bestk
## [1] 4
fit_asw$bestk
## [1] 4
plot(1:10,fit_ch$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio Calinski-Harabasz")

plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio silueta media")

Vamos ahora a realizar las gráficas de los modelos de 2, 3 y 4 clústeres para estudiar la agrupación de los datos en relación con el número de clústeres de forma visual.

Para el clusplot de dos clústeres, que es como tenemos agrupados los datos, podemos ver que aunque si existen dos grupos de datos claramente definidos, las elipses van en dirección contraria a la agrupación natural de los datos. Los mismo ocurre para los gráficos de 3 y 4 clústeres donde las elipses no definen claramente los núcleos de datos.

data.kmeans.2_clusters <- kmeans(x, 2)
clusplot(x, data.kmeans.2_clusters$cluster, color=TRUE, shade=TRUE, labels=2, lines=0)

data.kmeans.3_clusters <- kmeans(x, 3)
clusplot(x, data.kmeans.3_clusters$cluster, color=TRUE, shade=TRUE, labels=3, lines=0)

data.kmeans.4_clusters <- kmeans(x, 4)
clusplot(x, data.kmeans.4_clusters$cluster, color=TRUE, shade=TRUE, labels=3, lines=0)

Al contrastar los datos obtenidos en el análisis de reglas de semana del año y temperatura media podemos ver claramente la polarización de los datos, la cual se repite en los gráficos para 3 y 4 clústeres.

par(mfrow=c(2,2))
plot(data.kmeans.seasons[c(1,4)], col=data.kmeans.2_clusters$cluster)
plot(data.kmeans.seasons[c(1,4)], col=data.kmeans.3_clusters$cluster)
plot(data.kmeans.seasons[c(1,4)], col=data.kmeans.4_clusters$cluster)

El patrón se repite para la variable que define la vegetación de la zona, otra de las reglas relevantes obtenidas.

par(mfrow=c(2,2))
plot(data.kmeans.seasons[c(1,2)], col=data.kmeans.2_clusters$cluster)
plot(data.kmeans.seasons[c(1,2)], col=data.kmeans.3_clusters$cluster)
plot(data.kmeans.seasons[c(1,2)], col=data.kmeans.4_clusters$cluster)

data.kmeans.2_clusters.table <- table(data.kmeans.2_clusters$cluster,data.kmeans.seasons$outcome)
data.kmeans.2_clusters.table
##    
##     normal pandemic
##   1    615       98
##   2    421      322
data.kmeans.2_clusters.table.precision = 100*(max(data.kmeans.2_clusters.table[1], data.kmeans.2_clusters.table[3])+ max(data.kmeans.2_clusters.table[2], data.kmeans.2_clusters.table[4]))/(sum(data.kmeans.2_clusters.table))

data.kmeans.2_clusters.table.precision
## [1] 71.15385
data.kmeans.3_clusters.table <- table(data.kmeans.3_clusters$cluster,data.kmeans.seasons$outcome)
data.kmeans.3_clusters.table
##    
##     normal pandemic
##   1    272      301
##   2    476       98
##   3    288       21
data.kmeans.3_clusters.table.precision = 100*(max(data.kmeans.3_clusters.table[1], data.kmeans.3_clusters.table[4]) + max(data.kmeans.3_clusters.table[2], data.kmeans.3_clusters.table[5]) + max(data.kmeans.3_clusters.table[3], data.kmeans.3_clusters.table[6]))/(sum(data.kmeans.3_clusters.table))

data.kmeans.3_clusters.table.precision
## [1] 73.1456
data.kmeans.4_clusters.table <- table(data.kmeans.4_clusters$cluster,data.kmeans.seasons$outcome)
data.kmeans.4_clusters.table
##    
##     normal pandemic
##   1    278       21
##   2     62      124
##   3    227      178
##   4    469       97
data.kmeans.4_clusters.table.precision = 100*(max(data.kmeans.4_clusters.table[1], data.kmeans.4_clusters.table[5])+max(data.kmeans.4_clusters.table[2], data.kmeans.4_clusters.table[6])+max(data.kmeans.4_clusters.table[3], data.kmeans.4_clusters.table[7])+max(data.kmeans.4_clusters.table[4], data.kmeans.4_clusters.table[8]))/(sum(data.kmeans.4_clusters.table))

data.kmeans.4_clusters.table.precision
## [1] 75.41209

Al calcular la precisión de los tres modelos, aquel que hace uso de los 4 clústeres es el que mayor precisión presenta con un 75.41%. Este valor supera en un 4% y un 2% a los modelos de 2 y 3 clústeres respectivamente.

3.3.2 Uso de variables no estacionales

Repitamos ahora el procedimiento omitiendo la variable temporal de semana del año.

data.kmeans.no_seasons <- data.raw[, c(7,8,9,10,11,13)]
summary(data.kmeans.no_seasons)
##       ndvi        precipitation    air_temp_mean    diur_temp_range  
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:0.3239   1st Qu.:0.0255   1st Qu.:0.4005   1st Qu.:0.06621  
##  Median :0.4085   Median :0.0991   Median :0.5323   Median :0.10273  
##  Mean   :0.4356   Mean   :0.1172   Mean   :0.5375   Mean   :0.24174  
##  3rd Qu.:0.5353   3rd Qu.:0.1793   3rd Qu.:0.6864   3rd Qu.:0.42235  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000  
##     humidity          outcome    
##  Min.   :0.0000   normal  :1036  
##  1st Qu.:0.4397   pandemic: 420  
##  Median :0.6116                  
##  Mean   :0.5752                  
##  3rd Qu.:0.7157                  
##  Max.   :1.0000

Inicialmente observamos que la recomendación es de 2 clústeres.

x <- data.kmeans.no_seasons[,1:5]
d <- daisy(x) 

resultados <- rep(0, 7)
cluster  <- data.frame(c(2,3,4,5,6,7))

for (i in c(2,3,4,5,6,7))
{
  fit           <- kmeans(x, i)
  y_cluster     <- fit$cluster
  sk            <- silhouette(y_cluster, d)
  resultados[i] <- mean(sk[,3])
}

plot(2:7,resultados[2:7],type="o",col="blue",pch=0,xlab="Número de clusters",ylab="Silueta")

El método del elbow recomienda utilizar 3 clússteres para estos datos.

resultados <- rep(0, 7)
for (i in c(2,3,4,5,6,7))
{
  fit           <- kmeans(x, i)
  resultados[i] <- fit$tot.withinss
}
plot(2:7,resultados[2:7],type="o",col="blue",pch=0,xlab="Número de clusters",ylab="tot.tot.withinss")

Si procedemos con la función kmeansruns bajos los criterios de la silueta media y Calinski-Harabasz obtenemos que para ambos casos el número óptimo de clústeres es de 3 y 2 respectivamente, aunque para Calinski-Harabasz tenemos que los valores que validan la opción de 3 clústeres es similar a la de 2.

fit_ch  <- kmeansruns(x, krange = 1:10, criterion = "ch") 
fit_asw <- kmeansruns(x, krange = 1:10, criterion = "asw")
fit_ch$bestk
## [1] 3
fit_asw$bestk
## [1] 2
plot(1:10,fit_ch$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio Calinski-Harabasz")

plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio silueta media")

Procedemos pues a realizar el clusplot que claramente revela dos grupos y dos elipses, esta vez, en la dirección correcta.

data.kmeans.no_seasons.2_clusters <- kmeans(x, 2)
clusplot(x, data.kmeans.no_seasons.2_clusters$cluster, color=TRUE, shade=TRUE, labels=2, lines=0)

Al comparar la temperatura del aire media con otros factores, podemos ver una clara polarización de los datos cuando es contrastada con variables envueltas en las reglas de asociación, cosa que no ocurre con la precipitación.

par(mfrow=c(2,2))
plot(data.kmeans.no_seasons[c(1,3)], col=data.kmeans.no_seasons.2_clusters$cluster)
plot(data.kmeans.no_seasons[c(2,3)], col=data.kmeans.no_seasons.2_clusters$cluster)
plot(data.kmeans.no_seasons[c(4,3)], col=data.kmeans.no_seasons.2_clusters$cluster)

El model de tres clústeres se representa continuación.

data.kmeans.no_seasons.3_clusters <- kmeans(x, 3)
clusplot(x, data.kmeans.no_seasons.3_clusters$cluster, color=TRUE, shade=TRUE, labels=3, lines=0)

El contraste de variables se comporta similar de el dos clústeres.

par(mfrow=c(2,2))
plot(data.kmeans.no_seasons[c(1,3)], col=data.kmeans.no_seasons.3_clusters$cluster)
plot(data.kmeans.no_seasons[c(2,3)], col=data.kmeans.no_seasons.3_clusters$cluster)
plot(data.kmeans.no_seasons[c(4,3)], col=data.kmeans.no_seasons.3_clusters$cluster)

La precision es del 71.15%, pero en este caso para ambos modelos independienemete de los clústeres usados.

data.kmeans.no_seasons.2_clusters.table <- table(data.kmeans.no_seasons.2_clusters$cluster,data.kmeans.no_seasons$outcome)
data.kmeans.no_seasons.2_clusters.table
##    
##     normal pandemic
##   1    477       37
##   2    559      383
data.kmeans.no_seasons.2_clusters.table.precision = 100*(max(data.kmeans.no_seasons.2_clusters.table[1], data.kmeans.no_seasons.2_clusters.table[3])+ max(data.kmeans.no_seasons.2_clusters.table[2], data.kmeans.no_seasons.2_clusters.table[4]))/(sum(data.kmeans.no_seasons.2_clusters.table))

data.kmeans.no_seasons.2_clusters.table.precision
## [1] 71.15385
data.kmeans.no_seasons.3_clusters.table <- table(data.kmeans.no_seasons.3_clusters$cluster,data.kmeans.no_seasons$outcome)
data.kmeans.no_seasons.3_clusters.table
##    
##     normal pandemic
##   1    468       37
##   2    280      270
##   3    288      113
data.kmeans.no_seasons.3_clusters.table.precision = 100*(max(data.kmeans.no_seasons.3_clusters.table[1], data.kmeans.no_seasons.3_clusters.table[4])+ max(data.kmeans.no_seasons.3_clusters.table[2], data.kmeans.no_seasons.3_clusters.table[5]) + max(data.kmeans.no_seasons.3_clusters.table[3], data.kmeans.no_seasons.3_clusters.table[6]))/(sum(data.kmeans.no_seasons.3_clusters.table))

data.kmeans.no_seasons.3_clusters.table.precision
## [1] 71.15385

A pesar de que el modelo con los datos estacionales y cuatro clústeres denota un mejor rendimiento, consideramos que el de dos clústeres sin datos estacionales explica mejor el modelo. Esta conclusión la basamos en los siguientes hechos:

  • El modelo no estacional de dos clústeres se ajusta mejor a nuestro modelo de agrupación para los casos de dengue basado en la media.

  • El modelo no estacional hace uso de más variables para explicar las reglas mientras que el estacional delega la creación de reglas en la estación en la que se recogieron los datos.

  • Los datos estacionales son una constante independiente que sigue un patrón marcado y bien distribuido, cosa que no podemos afirmar del resto de las variables, que son medidas científicas de las condiciones climatológicas y morfológcas de las zonas estudiadas.

  • Los clusplot para los modelos estacionales invierte claramente la dirección de las zonas definidas y los clústeres de datos reales.

En resumen, creemos que la variable estacional explica el modelo desde una perspectiva poco científica y por ende será eliminada de los estudios que a continuación se producen.

3.4 Modelos no supervisados alternativos

Como método alternativo hemos seleccionado k-medoids basado en la distancia promedio de los puntos al punto real y central del cúmulo. Utilizaremos la función fviz_nbclust para realizar las estimaciones de los cúmulos. Calcularemos el número óptimos de clúster utilizando el método del codo y el de la distancia euclidiana. Como ya se ha mencionado anteriormente, omitiremos la variable estacional.

data.kmedoids <- data.kmeans.no_seasons
x <- data.kmedoids[,1:5]

El número recomendado de clústeres es de 2

fviz_nbclust(x = x, FUNcluster = pam, method = "silhouette", k.max = 10, diss = dist(x, method = "euclidean"))

fviz_nbclust(x = x, FUNcluster = pam, method = "wss", k.max = 10, diss = dist(x, method = "euclidean"))

El gráfico de clusplot es similar al del ejercicio anterior con dos grupos de datos claramente señalados.

data.kmedoids.cluster <- pam(x = x, k = 2, metric = "euclidean")
clusplot(x, data.kmedoids.cluster$cluster, color=TRUE, shade=TRUE, labels=2, lines=0)

La precisión del model es del 71.15%, igual que para el método de kmeans.

data.kmedoids.cluster.table <- table(data.kmedoids.cluster$cluster,data.kmedoids$outcome)
data.kmedoids.cluster.table
##    
##     normal pandemic
##   1    557      383
##   2    479       37
data.kmedoids.cluster.table.precision = 100*(max(data.kmedoids.cluster.table[1], data.kmedoids.cluster.table[3])+ max(data.kmedoids.cluster.table[2], data.kmedoids.cluster.table[4]))/(sum(data.kmedoids.cluster.table))

data.kmedoids.cluster.table.precision
## [1] 71.15385

A continuación vamos a utilizar el metido Fuzzy clustering que hace uso del algoritmo c-means donde el centro del clúster es la media de la probabilidad de las observaciones de pertenecer a ese grupo.

data.kmedoids <- data.kmeans.no_seasons
x <- data.kmedoids[,1:5]

Este método recomienda el uso de 3 clústeres.

data.fuz.output <- rep(0, 5)

for (i in c(2,3,4,5))
{
  fit           <- fanny(x = x, diss = FALSE, k = i, metric = "euclidean", stand = FALSE)
  data.fuz.output[i] <- fit$coeff['normalized']
}

plot(2:5,data.fuz.output[2:5],type="o",col="blue",pch=0,xlab="Número de clusters",ylab="Coeficiente Dunn")

data.fuz.3_cluster <- fanny(x = x, diss = FALSE, k = 3, metric = "euclidean", stand = FALSE)
clusplot(x, data.fuz.3_cluster$clustering, color=TRUE, shade=TRUE, labels=2, lines=0)

Podemos observar que el rendimiento es de 71.49%.

data.fuz.3_cluster.table <- table(data.fuz.3_cluster$cluster,data.kmeans.no_seasons$outcome)
data.fuz.3_cluster.table
##    
##     normal pandemic
##   1    302      122
##   2    256      261
##   3    478       37
data.fuz.3_cluster.table.precision = 100*(max(data.fuz.3_cluster.table[1], data.fuz.3_cluster.table[4])+max(data.fuz.3_cluster.table[2], data.fuz.3_cluster.table[5])+max(data.fuz.3_cluster.table[3], data.fuz.3_cluster.table[6]))/(sum(data.fuz.3_cluster.table))
data.fuz.3_cluster.table.precision
## [1] 71.49725

En general, concluimos que los modelos de dos clústeres observados devengan una precisión del 71.15% en todos los métodos, mientras que el de 3 clústeres para el método de Fuzzy clustering nos aporta un 71.49% de precisión, algo superior que el obtenido para kmeans.

3.5 Modelos supervisado “Decision Tree” sin PCA

Vamos a crear ahora un modelo basado en árbol de decisiones, para ello utilizaremos el paquete CARET que nos permitirá optimizar los parámetros del modelo de forma recursiva.

data.tree <- data.kmeans.no_seasons
summary(data.tree)
##       ndvi        precipitation    air_temp_mean    diur_temp_range  
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:0.3239   1st Qu.:0.0255   1st Qu.:0.4005   1st Qu.:0.06621  
##  Median :0.4085   Median :0.0991   Median :0.5323   Median :0.10273  
##  Mean   :0.4356   Mean   :0.1172   Mean   :0.5375   Mean   :0.24174  
##  3rd Qu.:0.5353   3rd Qu.:0.1793   3rd Qu.:0.6864   3rd Qu.:0.42235  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000  
##     humidity          outcome    
##  Min.   :0.0000   normal  :1036  
##  1st Qu.:0.4397   pandemic: 420  
##  Median :0.6116                  
##  Mean   :0.5752                  
##  3rd Qu.:0.7157                  
##  Max.   :1.0000

Antes de proceder, creamos nuestro conjunto de datos de entrenamiento y pruebas.

set.seed(777)
train <- createDataPartition(data.tree$outcome, p = 0.66, list = FALSE)

data.tree.train <- data.tree[train,]
data.tree.test <- data.tree[-train,]

data.tree.train.x <- data.tree.train[,1:5]
data.tree.train.y <- data.tree.train[,6]

data.tree.test.x <- data.tree.test[,1:5]
data.tree.test.y <- data.tree.test[,6]
data.tree.model.grid <- expand.grid(winnow = c(FALSE, TRUE), trials = seq(from = 1, to = 20, by = 2), model = c("rules"))

data.tree.model <- train(
  outcome ~ ., 
  data = data.tree.train, 
  method = "C5.0",
  tuneGrid = data.tree.model.grid,
  metric = 'Accuracy'
  )
data.tree.model
## C5.0 
## 
## 962 samples
##   5 predictor
##   2 classes: 'normal', 'pandemic' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 962, 962, 962, 962, 962, 962, ... 
## Resampling results across tuning parameters:
## 
##   winnow  trials  Accuracy   Kappa    
##   FALSE    1      0.7196039  0.2612001
##   FALSE    3      0.7128986  0.2730602
##   FALSE    5      0.7154400  0.2747109
##   FALSE    7      0.7152156  0.2705228
##   FALSE    9      0.7156527  0.2825051
##   FALSE   11      0.7170237  0.2830907
##   FALSE   13      0.7142482  0.2884227
##   FALSE   15      0.7151866  0.2892505
##   FALSE   17      0.7150693  0.2895645
##   FALSE   19      0.7148347  0.2886735
##    TRUE    1      0.7271019  0.2854781
##    TRUE    3      0.7185526  0.2882082
##    TRUE    5      0.7237425  0.2930462
##    TRUE    7      0.7252115  0.3020172
##    TRUE    9      0.7252247  0.3054888
##    TRUE   11      0.7246593  0.3034695
##    TRUE   13      0.7245092  0.3094721
##    TRUE   15      0.7251143  0.3091773
##    TRUE   17      0.7257747  0.3116507
##    TRUE   19      0.7254290  0.3092181
## 
## Tuning parameter 'model' was held constant at a value of rules
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were trials = 1, model = rules and winnow
##  = TRUE.

A continuación, extraemos el modelo más eficiente.

data.tree.model$bestTune
##    trials model winnow
## 11      1 rules   TRUE

Al siguiente gráfico nos revela la progresión de los modelos calculados según su precisión.

plot(data.tree.model)

Realizamos una predicción con los datos de pruebas, obsevando una precisión del 73.28%.

data.tree.predicted_model <- predict(data.tree.model, data.tree.test)
data.tree.predicted_model.confusionMatrix <- confusionMatrix(data.tree.predicted_model, data.tree.test$outcome)
data.tree.predicted_model.confusionMatrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction normal pandemic
##   normal      309       89
##   pandemic     43       53
##                                           
##                Accuracy : 0.7328          
##                  95% CI : (0.6914, 0.7713)
##     No Information Rate : 0.7126          
##     P-Value [Acc > NIR] : 0.1727          
##                                           
##                   Kappa : 0.2779          
##                                           
##  Mcnemar's Test P-Value : 8.975e-05       
##                                           
##             Sensitivity : 0.8778          
##             Specificity : 0.3732          
##          Pos Pred Value : 0.7764          
##          Neg Pred Value : 0.5521          
##              Prevalence : 0.7126          
##          Detection Rate : 0.6255          
##    Detection Prevalence : 0.8057          
##       Balanced Accuracy : 0.6255          
##                                           
##        'Positive' Class : normal          
## 

Podemos decir que este modelo realiza un mejor trabajo realizando predicciones del tipo normal en contraste con las pandemicas. Esto es un reflejo de la composición de los datos y su clasificación, donde existen más valores que dan como resultado un diagnóstico normal. Esto puedo verse reflejado en el alto nivel de sensibilidad ( true positive ) del test en comparación con la especificidad ( true negative ).

3.5.1 Modelos supervisado “Decision Tree” con PCA

A continuación, vamos a efectuar un análisis PCA para descartar variables que puedan estar aportando poco al modelo final. Este análisis nos permitirá simplificar el modelo con un bajo coste en la precisión final del modelo aunque con una ganancia en rendimiento a la hora de calcular el modelo a utilizar.

data.tree.pca.analysis <- prcomp(data.tree[,1:5], center = TRUE, scale. = TRUE)
summary(data.tree.pca.analysis)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5
## Standard deviation     1.3922 1.2904 0.8815 0.60091 0.50865
## Proportion of Variance 0.3876 0.3330 0.1554 0.07222 0.05174
## Cumulative Proportion  0.3876 0.7206 0.8760 0.94826 1.00000

Si observamos al resultado del modelo, podemos deducir que son los dos primeros componentes quienes explican mejor el modelo. De la rotación podemos observar que ndvi, air_temp_mean y diur_temp_range tienen más peso a la hora de calcular el modelo.

data.tree.pca.analysis$rotation
##                         PC1         PC2        PC3        PC4        PC5
## ndvi             0.61949453 -0.09773145  0.3302020 -0.5119503  0.4853335
## precipitation    0.23064993 -0.52108003 -0.7144128  0.2330397  0.3325403
## air_temp_mean   -0.41994894 -0.45286855  0.5264852  0.3566120  0.4628112
## diur_temp_range  0.61929755 -0.07972356  0.3026169  0.6311647 -0.3466537
## humidity        -0.05604084 -0.71237518  0.1087080 -0.3975655 -0.5652480

Si realizamos el plot de los dos primeros componentes, podemos observar que es PC1 quien explica mejor el modelo y que existe una mayor descentralización de los dos grupos. Además, podemos observar como ndvi_, air_temp_mean y diur_temp_range son las encargadas de polarizar el diagrama desde las esquinas éste.

ggbiplot(data.tree.pca.analysis,ellipse=TRUE, obs.scale = 1, var.scale = 1,  groups=data.tree$outcome, alpha = 0)

ggbiplot(data.tree.pca.analysis,ellipse=TRUE, choices = c(2,3),  obs.scale = 1, var.scale = 1,  groups=data.tree$outcome, alpha = 0)

Vamos ahora a repetir el modelo de predicción basado en el árbol de decisiones del paquete CARET.

data.tree.pca <- data.tree
summary(data.tree.pca)
##       ndvi        precipitation    air_temp_mean    diur_temp_range  
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:0.3239   1st Qu.:0.0255   1st Qu.:0.4005   1st Qu.:0.06621  
##  Median :0.4085   Median :0.0991   Median :0.5323   Median :0.10273  
##  Mean   :0.4356   Mean   :0.1172   Mean   :0.5375   Mean   :0.24174  
##  3rd Qu.:0.5353   3rd Qu.:0.1793   3rd Qu.:0.6864   3rd Qu.:0.42235  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000  
##     humidity          outcome    
##  Min.   :0.0000   normal  :1036  
##  1st Qu.:0.4397   pandemic: 420  
##  Median :0.6116                  
##  Mean   :0.5752                  
##  3rd Qu.:0.7157                  
##  Max.   :1.0000
set.seed(777)
train <- createDataPartition(data.tree.pca$outcome, p = 0.66, list = FALSE)

data.tree.pca.train <- data.tree.pca[train,]
data.tree.pca.test <- data.tree.pca[-train,]

data.tree.pca.train.x <- data.tree.pca.train[,1:5]
data.tree.pca.train.y <- data.tree.pca.train[,6]

data.tree.pca.test.x <- data.tree.pca.test[,1:5]
data.tree.pca.test.y <- data.tree.pca.test[,6]

Entrenamos un modelo con CARET indicando que queremos un preprocesamiento de los datos con PCA. CARET se encargará de seleccionar las variables óptimas para el modelo con variables reducidas.

data.tree.pca.model.grid <- expand.grid(winnow = c(FALSE, TRUE), trials = seq(from = 1, to = 15, by = 2), model = c("rules"))

data.tree.pca.model <- train(
  outcome ~ ., 
  data = data.tree.pca.train, 
  method = "C5.0",
  tuneGrid = data.tree.pca.model.grid,
  metric = 'Accuracy',
  preProcess = c('pca')
  )
data.tree.pca.model
## C5.0 
## 
## 962 samples
##   5 predictor
##   2 classes: 'normal', 'pandemic' 
## 
## Pre-processing: principal component signal extraction (5), centered (5),
##  scaled (5) 
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 962, 962, 962, 962, 962, 962, ... 
## Resampling results across tuning parameters:
## 
##   winnow  trials  Accuracy   Kappa    
##   FALSE    1      0.7067199  0.2181694
##   FALSE    3      0.7066854  0.1989197
##   FALSE    5      0.7092433  0.2144911
##   FALSE    7      0.7111387  0.2258201
##   FALSE    9      0.7105611  0.2289777
##   FALSE   11      0.7092170  0.2272623
##   FALSE   13      0.7094411  0.2278809
##   FALSE   15      0.7092170  0.2272623
##    TRUE    1      0.7098970  0.2082943
##    TRUE    3      0.7081524  0.2150357
##    TRUE    5      0.7102417  0.2178224
##    TRUE    7      0.7099807  0.2233506
##    TRUE    9      0.7115944  0.2348943
##    TRUE   11      0.7114823  0.2348682
##    TRUE   13      0.7117064  0.2354868
##    TRUE   15      0.7114823  0.2348682
## 
## Tuning parameter 'model' was held constant at a value of rules
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were trials = 13, model = rules and
##  winnow = TRUE.
data.tree.pca.model$bestTune
##    trials model winnow
## 15     13 rules   TRUE
plot(data.tree.pca.model)

Finalmente vemos que la precisión del modelo es de 71.46%, un 1.82% inferior que el anterior sin PCA.

data.tree.pca.predicted_model <- predict(data.tree.pca.model, data.tree.pca.test)
data.tree.pca.predicted_model.confusionMatrix <- confusionMatrix(data.tree.pca.predicted_model, data.tree.pca.test$outcome)
data.tree.pca.predicted_model.confusionMatrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction normal pandemic
##   normal      316      105
##   pandemic     36       37
##                                          
##                Accuracy : 0.7146         
##                  95% CI : (0.6725, 0.754)
##     No Information Rate : 0.7126         
##     P-Value [Acc > NIR] : 0.483          
##                                          
##                   Kappa : 0.1851         
##                                          
##  Mcnemar's Test P-Value : 1.024e-08      
##                                          
##             Sensitivity : 0.8977         
##             Specificity : 0.2606         
##          Pos Pred Value : 0.7506         
##          Neg Pred Value : 0.5068         
##              Prevalence : 0.7126         
##          Detection Rate : 0.6397         
##    Detection Prevalence : 0.8522         
##       Balanced Accuracy : 0.5791         
##                                          
##        'Positive' Class : normal         
## 

Este modelo ha aumentado su capacidad para realizar predicciones del tipo normal, pero ha perdido capacidad en lo referente a predicciones del tipo pandemic. Esto ha ocurrido como consecuencia de la eliminación de variables que, aunque aportaban poco al modelo, favorecían a la predicción en favor de la clase negativa o pandemic.

3.6 Conclusión árboles decisiones + PCA

Como bien se ha mencionado, existe un decremento en la precisión del modelo con PCA del 1.82%. Este es debido a que estamos trabajando con menos variables y por consecuente el modelo será menos preciso. Sin embargo, hay que tener claro que el propósito de PCA no es el de mejorar el rendimiento del modelo sino el de simplificarlo con el menor coste posible, que es lo que se ha logrado con este procedimiento.


4 Rúbrica